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Abstract 



We investigate tlie numerical stability of Cauchy evolution of linearized 
gravitational theory in a 3-dimensional bounded domain. Criteria of robust 
stability are proposed, developed into a testbed and used to study various 
evolution-boundary algorithms. We construct a standard explicit finite dif- 
ference code which solves the unconstrained linearized Einstein equations in 
the 3 -|- 1 formulation and measure its stability properties under Dirichlet, 
Neumann and Sommerfeld boundary conditions. We demonstrate the robust 
stability of a specific evolution-boundary algorithm under random constraint 
violating initial data and random boundary data. 
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I. INTRODUCTION 



The computational evolution of 3-dimensional general relativistic space-times by means 
of Cauchy evolution is a potentially powerful tool to study the gravitational radiation from 
black- hole/neutron-star binaries whose inspiral are expected to provide prominent signals to 
gravitational wave observatories. There are several 3-dimensional general relativistic codes 
under development to solve this problem. Boundary conditions are an essential part of 
these codes. At the outer boundary they must provide an outgoing radiation condition and 
extract the emitted waveform. For black-hole spacetimes, there is also an inner boundary, 
approximately given by the apparent horizon, where one excises the singular region inside 
a black hole. Instabilities or inaccuracies introduced at such boundaries have emerged as a 
major problem common to all code development. Historically, the first Cauchy codes were 
based upon the Arnowitt-Deser-Misner (ADM) formulation |l|J^ of the Einstein equations. 
Recently there has been pessimism that such codes might be inherently unstable because 
of the lack of manifest hyperbolicity in the underlying equations. In order to shed light 
on this issue, we present here a study of ADM evolution-boundary algorithms in the simple 
environment of linearized gravity, where nonlinear sources of physical or numerical instability 
are absent and computing time is reduced by a factor of five by use of a linearized code. 

Our two main results, for the case of fixed lapse and shift, are: 

• On analytic grounds, ADM boundary algorithms which supply values for all compo- 
nents of the metric (or extrinsic curvature) are inconsistent. 

• We present a boundary algorithm which allows free specification of the transverse- 
traceless components of the metric (or extrinsic curvature) at the boundary, and for 
which unconstrained, linearized ADM evolution can be carried out in a bounded do- 
main for thousands of crossing times with robust stability. 

The criteria for robust stability, which we present here, are the most severe that have 
been applied to Cauchy evolution in numerical relativity. The boundary algorithm differs 
from previous approaches and offers fresh hope for robust nonlinear ADM evolution. 

Our particular motivation for this work is the difficulty we have experienced imple- 
menting Cauchy-characteristic-matching (CCM) for 3-dimensional general relativity 
CCM provides a Cauchy boundary condition by matching the Cauchy evolution across the 
boundary to a characteristic evolution. For nonlinear scalar waves propagating in a flat 
3-dimensional space, CCM has been demonstrated to be more accurate and efficient than all 
other existing boundary conditions for Cauchy evolution 0, and it has been demonstrated 
mathematically that this conclusion also applies to gravity . In addition, in the spherically 
symmetric case of a self gravitating scalar wave satisfying the Einstein-Klein-Gordon equa- 
tions, CCM has been successfully applied at the inner boundary of a Cauchy evolution to 
excise the interior black hole region and, at the same time, at the outer boundary to provide 
a global evolution on a compactified grid extending to null infinity 0. These successes are 
promising for the application of CCM to 3-dimensional problems in general relativity but 
this has not yet been borne out. This difficulty, and the similar difficulty in efforts using 
perturbative matching P, may possibly arise from a pathology of the Cauchy boundary 
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which is independent of matching. In this work, we reveal such a pathology in the way 
boundary conditions have been applied in the ADM formulation of the Einstein equations 
which, at present, is the only formulation for which matching has been attempted. We also 
present a new form of ADM boundary algorithm which eliminates the pathology. 

The stability of the Cauchy evolution algorithm itself is straightforward to investigate 
by carrying out a boundary-free evolution on a 3-torus (equivalent to periodic boundary 
conditions). Such tests constitute Stage 1 of a 3-stage test bed for robust boundary stability 
which is summarized below and explained in detail in Sections |T|, |^ and 0. The periodic 
boundary tests serve to cull out algorithms whose boundary stability is doomed from the 
start. In earlier work, robust stability for characteristic evolution with random data on an 
inner boundary was demonstrated for characteristic evolution using the PITT Null Code p . 
In the course of the present investigation we have reconfirmed this robustness of the PITT 
code using the same specifications proposed here for Cauchy codes. 

CCM cannot work unless the Cauchy code, as well as the characteristic code, has a ro- 
bustly stable boundary. This is necessarily so because the interpolations between a Cartesian 
Cauchy grid and a spherical null grid continually introduce short wavelength noise into the 
neighborhood of the boundary. This is the rationale underlying the robustness criterion in 
our test bed. Robustness of the Cauchy boundary is a necessary (although not a sufficient) 
condition for the successful implementation of CCM. 

Analytic studies of Cauchy evolution of linearized gravity with boundaries at infinity 
reveal modes which grow linearly in time, but none which grow exponentially 0. The 
inaccuracy introduced by such secular modes can be controlled and is not of major concern, 
at least in the linearized theory. (Such secular modes can lead to exponential instabilities of 



numerical origin in the nonlinear theory if not properly treated [10|). In the case of a finite 



boundary, there is further potential for instability and a brief discussion is given in Sec. fT|. 

As is customary in numerical relativity, we monitor the existence of unstable modes by 
the growth of the Hamiltonian constraint. Because the constraints are not enforced during 
standard implementation of ADM evolution, the Hamiltonian constraint is an effective sensor 
of numerical instabilities. 

Stage 2 of the test bed is based on the simple boundary value problem obtained by 
opening one dimension of a 3-torus to form a 2-torus with plane boundaries normal to a 
Cartesian axis. Running a Cauchy-boundary algorithm with this topology and with random 
initial and random boundary data forms the second stage of our test bed, which is discussed 
in Sec. |V[ In Sec. [V], we present new evolution-boundary algorithms which are robustly 
stable. 

Stage 3 of the testbed is designed to test robustness of boundary conditions appropriate 
to an isolated system. In Sec. ^ we establish Stage 3 robustness of an ADM boundary 
algorithm. 

The main results presented here are experimental, in a computational sense. The difficul- 
ties encountered with finite Cauchy boundaries in general relativity have recently prompted 
some analytic investigations of the subject [|n|,0. However, these have so far been confined 
to hyperbolic formulations, as opposed to the ADM formulation, and to the analytic prob- 
lem, as opposed to the finite difference solution obtained by computation. Although it is 
not possible to make a direct comparison, the nature of our results are consistent with the 
general conclusions of these analytic studies. 
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There are several promising numerical approaches based upon hyperbolic (or "more 
hyperbolic") formulations of the equations [|T^ 23[|. Here we concentrate on ADM schemes, 



which are the most compact to code and require the least amount of memory because 
they have a smaller number of variables. Our results should provide useful benchmarks for 
other relativity codes. However, it should also be cautioned that the nature of a successful 
boundary algorithm is dependent on the form of the equations adopted, as well as the choice 
of discretization, and the ADM boundary algorithms we have obtained do not necessarily 
apply to other formulations. 

We use Greek letters for space-time indices and Latin letters for spatial indices. Four 
dimensional geometric quantities are explicitly indicated, such as ^"^^Rap and ^^^Gaf3 for the 
Ricci and Einstein tensors of the space-time, whereas Rij and R refer to the Ricci tensor 
and Ricci scalar of the Cauchy hypersurfaces. Linearized versions of these quantities are 
denoted by ^^^Raf3, Rij, etc. Three dimensional tensor indices are raised and lowered by the 
background Euclidean metric We write h = 6^^hij for 3-dimensional traces. We denote 
time derivatives by / = dtf. Our convention for the background Minkowski metric is such 
that the wave equation takes the form 

v"^dadp<!> = {-dl + dl + dl + dl)^ = 0. (1.1) 



II. GENERAL FRAMEWORK 

A. The linearized ADM system 

The ADM formulation of the Einstein equations introduces a foliation of space-time by 
a time coordinate t and expresses the 4-dimensional metric as 

ds^ = -a^dt^ + gij {dx' + (3'dt) {dx^ + (3^dt) , (2.1) 

where Qij is the induced 3-metric of the t = const slices, a is the lapse and /3* the shift, 
with the normal to the slices given by = (1, —(3^) /a. The equations '"^^Rij = yield the 
evolution equations 

dtQij - £pgij = -2aKij (2.2) 
dtKi, - £fsKy = -DiD.a + a {Rij + KKij - 2K\Kij) , (2.3) 

for the 3-metric Qij and the extrinsic curvature Kij = —^£ngij, subject to the constraints 

R - KijK'^ + K^ = (2.4) 
Dj{K'^-g'^K) =0. (2.5) 

Here R, Rij and Di are the Ricci scalar, Ricci tensor and connection of the 3-metric, respec- 
tively. 

For simplicity we consider a gauge in which the lapse is unity and the shift vanishes 
(Gaussian coordinates), so that the linearized metric gais = rjap + hap satisfies hta = 0, and 
obeys the linearized ADM evolution equations 
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dthij — ~2Kij 

dtK,, = Ri,, (2.6) 
subject to the (linearized) constraints 

^ = 

dj{K'^-6'^K) =0. (2.7) 

Here we consider a 1-parameter system of equations, equivalent to the linearized Einstein 
equation, consisting of the six evolution equations Sij = along with the four constraint 
equations C = Cj = 0, where 

S,, := + hSif, (2.8) 

C := ^^^Gtt, Ci := —^^^Gti and the parameter A allows mixing the (linearized) Hamiltonian 
constraint C into the evolution equations. For A = we recover the standard ADM system. 

Codes under development for the evolution of 3-dimensional space-times without sym- 
metry apply the constraint equations at the initial time but do not enforce them during the 
evolution. It is crucial for this approach that the constraints be stably propagated in time. 
An investigation by Frittelli shows that this requires the parameter A in Eq. (|2.8|) satisfy 



1 + A > 0. This follows from an analysis of the linearized Bianchi identities d^^^^G^ = 0, 
which imply that 

C + + A) d'C + djS'^ = (2.9) 
C + diC = 0. (2.10) 

Thus if the evolution equations are satisfied then the Hamiltonian constraint satisfies 

C - (1 + A)9^9fcC = 0; (2.11) 

This equation has a well-posed initial value problem for A > — 1 (when it is hyperbolic) and 
also for A = — 1, but for A < — 1 the equation is elliptic and the initial value problem is not 
well-posed. In the standard ADM case, A = and the Hamiltonian constraint propagates 
along the light cone. We consider here evolution equations with a range of A. 
The linearized evolution equations ( p2.8| ) take the form 

2Kij 

Ki, = -\d^d^K, + ^{diH, + d,Hi) + ^5,,AC, (2.12) 

where 

H, = d^{h,-U,jh). (2.13) 
and we can express the Hamiltonian as 

C = ]^d,W--^dmd'^h. (2.14) 
A spectral analysis of a system similar to Eq's ( p.l2[ ) - ( |2.13| ) is presented in . 
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B. Finite difference algorithms 



The evolution variables consist of the 3-metric perturbations hy and their associated 
momentum Kij = —hij/2. The evolution is implemented on a uniform spatial grid 
{xj,yk,zi) = {jAx,kAx,lAx) with time levels = nAt. The three different evolution 
algorithms we apply can be discussed in reference to the scalar wave Eq. ( |1 . 1| ) , rewritten in 
the form 

$ = -2K 

k=~dmd"'^, (2.15) 

analogous to the first differential order in time and second differential order in space form 
of the ADM equations. We denote ^^f^i = jAx, fcAx, /Ax). All second derivatives on 
the right hand side of Eq. ( |2.15| ) are calculated as centered 3-point finite differences. 



1. Standard leapfrog (LFl) 

The first evolution algorithm, which we refer to as LFl, is a standard leapfrog imple- 
mentation of Eq. ( p. 15 ): 



where is the second order accurate centered difference approximation to the Laplacian. 
It is known that this algorithm has a time-splitting instability in the presence of dissipative 



and nonlinear effects 26 



2. Staggered leapfrog (L'F2) 

The second evolution algorithm, which we refer to as LF2, is a staggered in time leapfrog 
scheme which is not subject to the time-splitting instability: 

1 



^Tk:i = nk,i-'^K-jAt (2.17) 



■j,k,l — ^j,k,l "^^j^k^l 

KtkJ = K^uJ -l^'n^A^- (2-18) 
Here K is evaluated on the half grid. Subtraction of the equation 



from Eq. ( p.l7| ) and elimination of K using Eq. ( p.l8| ) shows that LF2 is equivalent to the 



standard centered second-order scheme for the second differential order in time form of the 
wave equation ( |1 . 1| ) , in which $ lies on integral time levels and K is not introduced. 
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3. Iterative Crank-Nicholson (ICNj 



The third evolution algorithm, which we refer to as ICN, is a two-iteration Crank- 
Nicholson algorithm. For an A^-iteration Crank-Nicholson algorithm, the following sequence 
of operations is executed at each time-step: 

1. Compute the first order accurate quantities 



(0).,i 1 
2 



Ktk] = Kl,,i - ^V'^kA^- (2-20) 



2. Starting with z = 0, compute the midlevel values 

AmI = a'SI ^ (2-21) 



3. Update using levels n and n + |, 

A-'»tl = A'".,, - l'-|Ai. (2.22) 

4. Increment i by one and return to step ^ until z = is reached. 

A stability analysis by Teukolsky shows that the evolution scheme is stable for A^ = 2 
and A^ = 3 iterations, unstable for A^ = 4 and N = 5, stable for A^ = 6 and N = 7, etc. IE? . 



4- The Courant-Friedrichs-Lewy condition 



The stability of the three evolution algorithms requires they obey the Courant-Friedrichs- 
Lewy (CFL) condition that the numerical domain of dependence contain the analytic domain 
of dependence, a common requirement for explicit algorithms. For the staggered leapfrog 
(LF2) and the iterative Crank-Nicholson (ICN) algorithms, we set At = Ax/4 and for the 
standard leapfrog (LFl) we set At = Ax/8, in all cases slightly less than half the CFL 
condition for the algorithm. 
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5. Boundary conditions 



Boundary conditions are implemented computationally in the following way, which we 
illustrate in terms of scalar wave boundary data specified at 2; = in terms of a function /. 
The Dirichlet condition 

$(t,x,i/,0) = /(t,x,i/) (2.23) 
is straightforward to implement as 

^KO = fh- (2.24) 

The Neumann condition 

a,$(t,x,i/,0) = /(t,x,y) (2.25) 
is implemented as a 3-point one-sided derivative 

*i,^,o = ^(-"^",^,2 + ^^Ki - '^^^fW- (2-26) 
The Sommerfeld condition 

{dt-d,)^{t,x,y,^) = f{t,x,y), (2.27) 



is implemented in the interpolative form used in several relativity codes |T3| , pr5| , pT| by mod- 
eling the field in the neighborhood of the boundary as $(t + z,x,y) and using a 3-point 
spatial interpolation to obtain 

- U2 - ^ Vl - — ^ + — ^2 - 

^jAo - 2 Axj V ^'''^ Ax Ax) ^'^'1 



III. STAGE 1: ROBUST EVOLUTION STABILITY 

Periodic boundary conditions are equivalent to a toroidal topology and do not introduce 
the local effects of a real boundary. They provide a test of the evolution code isolated from 
the effects of boundary conditions. Because an instability in such a code may not be evident 
for a considerable time if masked by a strong initial signal, the use of random data is efficient 
at revealing instabilities early in the evolution. Random initial data does not satisfy the 
constraints but that poses no difficulty here, where we are only concerned with stability. 
These observations motivate our Stage 1 test bed: 

Stage 1: Run the evolution code on a 3-torus with random initial Cauchy data. 
The stage is passed if the Hamiltonian constraint C does not exhibit exponen- 
tial growth. 
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An evolution code which does not exhibit exponential growth under these conditions 
is defined to be robustly stable. Failure at Stage 1 would rule out applications with 
boundaries. 

We use an evolution time of 2000 crossing times (2000L, where L is the linear size of the 
computational domain) on a uniform 48^ spatial grid with a time step slightly less than half 
the Courant-Friedrichs-Lewy limit. These conditions are computationally practical and are 
used to determine whether there is exponential growth of the Hamiltonian, as measured by 
the £oo norm. All runs reported in this paper are made with these specifications. 



A. Stage 1 results 

We have applied the Stage 1 test to determine whether any of three evolution algorithms, 
LFl, LF2 and ICN are intrinsically unstable. We apply the test on the fiat 3-torus de- 
termined by the periodicity conditions hij{x,y,z) = hij{x + L,y,z) = hij{x,y + L,z) = 
hij{x,y,z + L). The Cauchy data, {hij,hij), can be initialized as random numbers in any 
interval {—A, A), since the system is linear. Here we use the interval (— 10^^, +10~^). 

When applied to the scalar wave equation in the hybrid first order in time and second 
order in space form of Eq. ([2.15|), all three algorithms LFl, LF2 and ICN pass Stage 1. 



This confirms, for the case of random data, prior work showing that the hybrid system 
has a well behaved computational evolution. 



Furthermore, when applied to the ADM system ( |2.12| ) for gravitational evolution with 
A equal to 0, 2 and 4, all three algorithms LFl, LF2 and ICN also pass Stage 1. For 
runs with A equal to -0.1, 4.1 and 5.0 these three evolution algorithms exhibited exponential 
growth. This indicates a range of stability for < A < 4. In this range, it is notable that 
the norm of the Hamiltonian constraint grows linearly in time for LFl and LF2 but decays 
exponentially for ICN. This apparently results from the artificial dissipation inherent in 
ICN. 

The stability of the discretized system of ADM equations is more restrictive than (but 
consistent with) the range — 1 < A found by Frittelli for stable evolution of the con- 



straints in the continuum theory. For A = —1, algorithms LFl and LF2 show exponential 
growth whereas the norm of the Hamiltonian only grows linearly for ICN. However, for 
A = —1.01 or A = —0.99 this norm grows exponentially for ICN. This anomalous behavior 
suggests that the the special case A = — 1 (the Einstein system of evolution equations) can 
be successfully evolved but that its numerical stability is highly sensitive to the choice of 
finite difference scheme. For that reason, we have not investigated this case in the presence 
of a boundary. The upper limit of the window of stability at A = 4 is related to the size of 
the time step. For algorithm LF2, a run with At = Ax/8 (half the time step of the standard 
runs) and A = 20 showed no exponential growth. This seems to arise from the increase of 
the constraint propagation speed with A, which makes the Courant-Friedrich-Lewy condition 
more stringent. 

In summary, the hybrid scalar wave system ( |2.15|) and the ADM system (|2.12| ) , 



with A equal to 0, 2 and 4, pass Stage 1 for the three evolution algorithms LFl, 

LF2 and ICN. These are the evolution systems whose boundary stability we investigate in 



Sec. IV 
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IV. STAGE 2: ROBUST BOUNDARY STABILITY 



The general linear hyperbolic equation in second order differential form for a scalar 
field has a well posed Cauchy problem in a region with Dirichlet, Neumann or Sommerfeld 
boundary conditions (e.g., see p9|). For a system of coupled scalar fields, or a tensor 
field with coupled components, it is standard practice to reduce the equations to first order 



differential form in order to examine hyperbolicity and appropriate boundary conditions 
For a first order system in diagonalizable, strongly hyperbolic form there is a straightforward 
way to decide which variables require data at a given boundary |^ . Variables propagating 



along future directed characteristics which emanate from the boundary can be assigned 
free data, but assigning arbitrary boundary values to the remaining variables would be 
inconsistent with the evolution equations. When a second order system is reduced to first 
order form, spatial derivatives of the field become auxiliary variables, so that there is no 
longer any natural distinction between, say, Dirichlet and Neumann boundary conditions. 
What might have been termed a Neumann condition in the original system now appears in 
Dirichlet form. Friedrich and Nagy |]12| have recently given a complete treatment of a well- 
posed boundary-initial-value problem for a symmetric hyperbolic version of the nonlinear 
vacuum gravitational equations. They find a continuum of allowed boundary conditions on 
the field variables, which can be more naturally distinguished as ranging from "electric" to 
"magnetic" . 

The hybrid form of the scalar wave equation( |2.15| ) does not not fit into any hyperbolic 
category but, since the spatial derivatives of the field are not treated as auxiliary variables, 
we retain the classification of Dirichlet, Neumann and Sommerfeld boundary conditions. 
Following common practice, we also retain this classification in the case of the ADM system 



Whereas Stage 1 tests stability of the interior evolution algorithm itself. Stage 2 is de- 
signed to be a simple stability test of the combined evolution-boundary algorithm. The 
boundary algorithm by itself is neither stable nor unstable; rather the combination of the 
boundary algorithm with a (stable) evolution algorithm may be stable, and a combination 
with another (stable) evolution algorithms may be unstable P^. In Stage 2, the three torus 



is opened up in the z-direction to form a space of topology x [0, L], with boundaries at 
z = and z = L coinciding with planes of grid points. A boundary algorithm for these 
points is necessary in order to update the evolution at grid points neighboring the boundary. 
One purpose of the testbed is to measure suitability for matching the Cauchy evolution to 
an exterior numerically generated solution, such as in CCM, where interpolations between 
the exterior and interior grids continually introduce random error at the Cauchy boundary. 
This motivates the following criterion for robust evolution-boundary stability: 

Stage 2: Run the evolution-boundary code on x [0, L\ with random initial 
Cauchy data and random boundary data. The stage is passed if the Hamilto- 
nian constraint C does not exhibit exponential growth. 

As an illustration of how Stage 2 is implemented, rather than giving smooth Dirichlet 
data, such as the homogeneous data $(t,a;, |/,0) = for a scalar field, we require that $ 
be prescribed as a random number at each boundary point. Similarly, in the Neumann or 
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Sommerfeld cases, (9^$ or {dt — d^)^ are prescribed as random numbers at z = 0. In order 
to avoid inconsistencies, the initial and boundary data are both set to in a few grid zones 
near the intersection of the initial Cauchy surface with the boundary. 

As a first set of experiments, we have confirmed that the hybrid scalar system (p.l5|) 
passes Stage 2 for the three evolution algorithms LFl, LF2and ICN with a Dirichlet 
boundary algorithm. Sommerfeld and Neumann boundary algorithms were less success- 
ful, as indicated in Table 1. Only those combinations of evolution and boundary algorithms 
which are robustly stable for a scalar field should be expected to pass Stage 2 for the ADM 
system. 

Next, we tested ADM evolution with boundary data prescribed for each component of 
Kij, as has been common practice. In Sec. 0, for the case A = 0, we show this practice leads 
to an inconsistent evolution-boundary problem, whose finite difference solutions cannot in 
general converge to a correct continuum solution. It is notable that the numerical results 
were quite mixed, not necessarily showing unstable growth. For homogeneous boundary 
conditions and evolution with A = 0, we first found that the only stable combination was 
ICN evolution with a homogeneous Sommerfeld boundary. (All other combinations showed 
exponential growth on the order of 10 crossing times.) Next, we applied random Sommerfeld 
boundary data to all components of Kij (the analogue of choosing / randomly in Eq. ( |2.27D ), 
again with ICN evolution and A = 0. The log plot in Fig. |l| shows the Hamiltonian constraint 
growing at late times as t", for n ^ 1.92. Such polynomial growth is normally regarded as 
stable. However, in this case, there is a large multiplying constant, and the magnitude of 
the error (of the order of 1000 at t = 2000) is unacceptably high. 



V. NEW ADM BOUNDARY ALGORITHMS 

A. Consistency of ADM boundary conditions 

Various types of boundary conditions can be applied to a scalar wave, e.g. Dirichlet, 
Neumann or Sommerfeld. There are more options in the ADM case, corresponding to 
Dirichlet, Neumann or Sommerfeld conditions on the various components of the metric, 
or equivalently on the components of extrinsic curvature Kij. However, many of these 
versions are inconsistent with the evolution or constraint equations. In this regard, we list 
some combinations of the linearized Einstein equations and their implications for a correct 
boundary algorithm. We take the boundary to be a surface z = const and denote the 
transverse directions by = {x,y). 

• The linearized Einstein equation component 

2 (^)G^ = 2ki + d^dBhi - dAdBh"^^ = 0, (5.1) 

can be apphed on the boundary to evolve the transverse trace = K^^ + Kyy, given 
the transverse-tracefree components Kab — ^SabKq. 

• The linearized Ricci tensor equation 
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(5.2) 



can be applied on the boundary to evolve the trace K, thus determining K^z in terms 
of transverse components. 

• The Einstein equation components 

2 = -2k^ - dnid'^hf - a^/if ) - dzd^'h^ + dzd^'hi = 0. (5.3) 

can be applied on the boundary to determine K^, given the transverse components. 

• The linearized momentum constraint 

= dzK^' + SbK^^ - d^K = (5.4) 

or the combination of the time derivative of the momentum constraints with Eq. ( ^.21) , 

give other ways to update the Neumann boundary data for dzK^ in terms of Dirichlet 
boundary values of Kab- 

• The combination 

- d'^^^R\ = dzkz, + OaK^ = (5.6) 

can be used to update the Neumann boundary data for K^z- 

In the symmetric hyperbolic treatment of the Einstein equations by Friedrich and 
Nagy ||I2|, only 2 components of the Weyl tensor can be prescribed as free boundary data. 
It would thus be surprising if free boundary values could be assigned to all metric variables 
(or their associated momentum variables) for an ADM system with gauge freedom fixed by 
an explicit choice of lapse and shift, as shown by the following proposition. 

Proposition: Prescription of Dirichlet boundary data on all components of the metric 
(or extrinsic curvature) of the ADM system (|2.12|) with A = gives rise to an inconsistent 
evolution-boundary problem. The same is true for Neumann or Sommerfeld boundary data. 

Proof: Consider homogeneous Dirichlet data consisting of setting all components of hij to 
zero on the boundary. Then the function \1/ := dAd^h^ — d'^d^hAB vanishes on the boundary 
and Eq. ( |5.3| ) (one of the evolution equations for this system) implies that the normal 
derivative dz"^ also vanishes on the boundary. But it is easy to verify, in the case A = 0, 
that the evolution equations for hij imply that \l/ satisfies the scalar wave equation. Thus 
the vanishing Dirichlet data for hij generates, for any initial data, a solution of the wave 
equation whose Dirichlet and Neumann boundary data both vanish. This a classic example 
of an inconsistent boundary value problem for the scalar wave \I'. (This is evident from 
considering the fate of an initial pulse of compact support when it reaches the boundary.) 



12 



Similarly, consider the homogeneous Sommerfeld data {dt — dz)hij = applied to a// met- 
ric components on the plane boundary at 2; = 0. If hij were a global solution consistent with 
this boundary data then, since the equations are linear and have space-time translational 
symmetry, hij = {dt — dz)hij would also be a global solution but with vanishing Dirichlet 
data for all components at the boundary. Thus, as in the Dirichlet problem, a Sommerfeld 
boundary condition, or by the same argument a Neumann boundary condition, applied to 
all components of the metric also leads to an inconsistent boundary value problem. □ 



B. Robustly stable Dirichlet evolution-boundary algorithms 

In order to formulate consistent boundary algorithms, we denote by hrpT the traceless 
part of the components transverse to the boundary, i.e. (/i^^; — hyy) and /i^^ in our Stage 
2 test with boundaries at z = and z = L. Since our gauge choice /if^ = is consistent 
with the radiation gauge subclass of harmonic coordinates, these TT components represent 
the free modes of waves propagating normal to the boundary. We make the hypothesis 
that the boundary values of hxr, or equivalently K^t, should be freely specified in either 
Dirichlet, Neumann or Sommerfeld form. This is motivated in the Dirichlet case by the 
consistency of characteristic evolution where the free data on a worldtube corresponds to 
Dirichlet data for Htt in the linearized approximation. Given this TT boundary data, the 
boundary algorithm must determine boundary values of the remaining components using 
the linearized gravitational equations. 

The following five Dirichlet boundary algorithms exhibit Stage 2 robust sta- 
bility for the ICN evolution algorithm. The algorithms update the boundary values 
of the extrinsic curvature, with boundary values for the metric perturbation updated by 
the centered difference version of the first of Eq's ( ^.12[ ). Given random initial and bound- 



ary data for the transverse-traceless components Ktt, ah five boundary algorithms update 
the boundary values of the trace via integration of Eq. ([Ol) . Boundary values of the 
remaining unspecified components are updated as follows: 



BAl: We apply Eq. ( p.6| ) to update K^z and Eq. ( p. 3D to update K^. 



BA2: We apply Eq. (|5.2|) to update Kzz and the momentum constraint Eq. ( ^.4|) 
to supply boundary values for dzK^ which, expressed as a 3-point sideways finite 
difference, are used to update . 



BA3: We apply Eq. ( |5.2|) to update boundary values for Kzz and Eq. ( ^.5|) to supply 
boundary values for dzKf which, as in BA2, are used to update using a centered 
time difference. 



BA4: We apply Eq. (|]2D to update Kzz and Eq. (|]3D to update K^. 



BA5: We apply Eq. ( |5l6| ) to update Kzz (with finite difference stencils as above) and 
Eq. (^) to update K^. 



All five boundary algorithms satisfy Stage 2 robust stability for ICN evolution. Fig. |^ 
shows the behavior of the Hamiltonian constraint for these five algorithms in the case A = 2. 
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Note that BA2 and BA3 have identical performance, as might be expected as they differ 
only with respect to details of initialization at the boundary. BAl and BA4 show less noise 
in the Hamiltonian constraint than the others, with BA5 showing the largest (although still 
linear) growth. BAl gave the best performance, with the Hamiltonian constraint actually 
decreasing slowly at late times. 

For A = and A = 4, boundary algorithms BAl and BA2 are also robustly stable for 
ICN evolution. (The other boundary algorithms were not checked for these cases in order 
to conserve computing time). 

While these 5 Dirichlet boundary algorithms were robust for ICN evolution, they failed 
Stage 2 for LFl and LF2 evolution with A = 0, 2 and 4, with the exponential growth 
rate typically decreasing with increasing A. Table 2 summarizes the performance for A = 2. 
The failure of these evolution-boundary algorithms for leapfrog evolution, but not for ICN, 
emphasizes the complexity of the finite difference problem compared to the corresponding 
analytic problem. 

C. Neumann and Sommerfeld boundaries 

We attempted to modify the Dirichlet boundary algorithms BAl to BA5 to obtain 
stable evolution with Neumann or Sommerfeld boundary data specified for the extrinsic 
curvature components Ktt- In the Neumann case, assuming all components of the metric 
have been determined at time level — 1 and the evolution has been apphed to update 
all components at level N except at the boundary, we express the Neumann boundary data 
dzK^T in finite difference form according to Eq. ( |2.2(j| ) to update Ktt on the boundary at 
level A^. This supplies the necessary data to apply the Dirichlet boundary algorithms to 
update all remaining components. 

Similarly, in the Sommerfeld case, given that the metric has been determined at level 
— 1 and the evolution has been applied to update all components at level N except at 
the boundary, we apply the interpolative Sommerfeld condition in finite difference form 
according to Eq. (|2.28| ) to update Ktt on the boundary at level N . Again this supplies 
the necessary data to apply the Dirichlet boundary algorithms to update all remaining 
components. 

We tried an extensive, although not exhaustive, set of combinations of evolution al- 
gorithms, boundary algorithms and values of A with Sommerfeld or Neumann conditions 
applied to the TT components but we were unable to obtain acceptable Stage 2 evolution. 

VI. TESTS WITH A CUBIC BOUNDARY (STAGE 3) 

For application of these algorithms to an isolated astrophysical system, we next perform 
tests with a cubic boundary. This is the standard boundary geometry adopted in Cauchy 
evolution codes based upon Cartesian coordinates. We propose the following operational 
criteria of robust stability for a Cartesian evolution-boundary algorithm for an isolated 
system: 
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Stage 3: Run the evolution-boundary code with a cubic boundary with random 
initial Cauchy data and random boundary data. The stage is passed if the 
Hamiltonian constraint C does not exhibit exponential growth. 

In view of the Stage 2 results, we confine our Stage 3 investigation to ICN evolution with 
Dirichlet boundary data on all faces of the cube applied with the (best performing) boundary 
algorithm BAl . The edges and corners of the cube must be handled separately. The two 
components K^t = —\hTT are treated as free data (i.e. are specified randomly) on all faces, 
edges and corners. While this means two free quantities and four update equations on the 
faces, there are four free quantities on the edges so that one only needs two update equations. 
Furthermore, at any corner, the six Ktt components from the neighboring faces, K^y, K^z, 
Kyz, Kxx — Kyy, K^x " K^z 8.11(1 Kzz " K yy , oxQ reduccd to five that are independent and 
therefore freely specifiable by means of the identity [Kxx—Kyy] + [Kyy—Kzz] + [Kzz—Kxx] = 0. 
Thus only one equation is needed to update the corners. 

As just indicated, all non-diagonal components are freely specified TT data on the cor- 
ners. Given the additional TT data [K^x — Kyy] and [fC^^ — K^x], the missing diagonal 
component Kxx is computed from 

S-ft'xx = K -\- [Kxx Kyy] + [Kxx Kzz] ) 

where K is updated using the equation 

= -k = o. 

It remains to give the algorithm for the edges. On the edges parallel to the x-axes, 
Kxy and Kxz are specified as boundary data. The missing non-diagonal component Kyz is 
updated using ^^^Gyz = 0, the same equation used on the neighboring faces except now the 
derivatives of the metric in the y and z directions must be computed by 3-point sideways 
differencing. The diagonal components of Kij are computed the same way as on the corners. 

We should note that the routine that solves the constraint 

(-C" + 9"(^)i?J) = (6.1) 

on a face of the cube with normal in the n-direction must be called after the missing non- 
diagonal components have been updated on the edges surrounding that face. Otherwise, in 
the case of the z = const face, when computing the quantity dyKyz on the top time-level, 
with centered finite differencing, one might use values of Kyz on the edge parallel to the 
X-axis that were not yet updated. 

We confirmed that the above algorithm is robustly stable by performing runs with 
A = 0,2,4 and random initial and boundary data. The behavior of the Hamiltonian con- 
straint function of time is shown in Fig. H. 

VII. CONCLUSION 

We have shown that linearized ADM evolution with boundaries can be carried out with 
long term stability in a test bed consisting of random constraint violating initial data and 
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random boundary data applied to the trace-free-transverse metric. Adding the Hamiltonian 
constraint (with A > 0) to the Ricci system of hnearized equations appears to give better 
performance, but does not drastically affect overall robustness. The successful implementa- 
tion of an ADM boundary algorithm presented here offers new hope both for the long term 
stability of nonlinear ADM evolution and for the prospects of matching an exterior solu- 
tion at an ADM boundary. However, this optimism should be tempered with the following 
caveats. 

Although we have tested our algorithm only for the case of unit lapse and zero shift, 
the extension to any explicitly assigned values of the lapse and shift appears to be straight- 
forward, at least in the linearized theory. However, the use of dynamical gauge conditions, 
which couple the values of lapse and shift to the metric, would require case-by-case recon- 
sideration. 

A spherical boundary would be necessary for an application of our algorithms to CCM, 
The implementation of a spherical boundary algorithm is simple in principle. Only the 
TT metric (or extrinsic curvature) components should be matched at the boundary, with 
the the remaining components updated using the evolution equations. The identification 
of the TT components can be readily made in the local tangent space of the boundary. 
However, a preliminary investigation reveals nontrivial technical problems arising from the 
non-alignment of a spherical boundary with the Cartesian grid. ||33|| . 

Results for the linear theory are important for ruling out approaches that cannot work in 
the nonlinear case. However, the real value of our robust boundary algorithms will depend 
upon whether they can successfully be applied to the nonlinear ADM equations. For BAl, 
the best performing boundary algorithm, the formal implementation appears to be straight- 
forward (when the lapse and shift are given explicitly). It is standard practice in numerical 
evolution to choose the boundary to follow the evolution, so that the grid is propagated up 
the boundary. This allows straightforward identification of the TT components of the metric 
and extrinsic curvature. An examination of the nonlinear equations used in the boundary 
algorithm shows that second derivatives do not appear in any essentially new way that would 
alter the finite difference stencils. Nonlinear terms with first time derivatives which appear 
in the boundary update scheme can be evaluated either by means of iterative techniques or 
in terms of previously known time levels by backwards differencing. A separate and more 
problematic issue is the stability of such an implementation. At the very least, stability in 
the nonlinear case would require suppressing the secular modes of the linear theory from 
becoming exponential |TD[. Preliminary work underway to incorporate our boundary 



algorithms in a nonlinear ADM code shows improved performance in the weak field regime 
over applying boundary conditions to all components of the metric, but it is premature to 
judge robust nonlinear stability. Our results for the linearized equations could not have been 
obtained without substantial computational experimentation and the same certainly holds 
for their extension to the nonlinear case. 
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TABLES 



TABLE I. Results of Stage 2 tests for scalar wave evolutions with the LFl, LF2 and ICN 

algorithms, using Dirichlet, Sommcrfeld and Neumann boundary conditions on the z = const faces 
of a cube. A "-y indicates robust stability, a "x" indicates exponential instability and a "?" 
indicates non-linear growth which was not clearly exponential on the time scale of the test. 





Dirichlet 


Sommerfeld 


Neumann 


LFl 




X 


? 


LF2 






? 


ICN 






? 



TABLE IL Stage 2 tests of ADM evolution with A = 2 for boundary algorithms BAl - BA5. 
A "-y/" indicates robust stability. A "x" indicates instability with the exponential growth rate 
indicated in units of crossing time (CT). 





BAl 


BA2 


BA3 


BA4 


BA5 


LFl 


X (150 CT) 


X (300 CT) 


X (300 CT) 


X (300 CT) 


X (300 CT) 


LF2 


X (25 CT) 


X (300 CT) 


X (100 CT) 


X (25 CT) 


X (100 CT) 


ICN 
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FIGURES 
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time / L 

FIG. 1. A log plot of the Iqo norm of the Hamiltonian constraint as a function of crossing time 
for a Stage 2 test of random Sommerfeld boundary conditions on all metric components. 
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FIG. 2. Stage two performance of the Hamiltonian constraint as a function of crossing time for 
the five robustly stable algorithms BAl to BA5 . 
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FIG. 3. Behavior of the Hamiltonian constraint for a Stage 3 test with cubic boundary. 
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